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Abstract 

We completed the development of simulation code that is designed to study 
the behavior of a conjectured dark matter galactic halo that is in the form 
of a Bose-Einstein Condensate (BEC). The BEC is described by the Gross- 
Pitaevskii equation, which can be solved numerically using the Crank-Nicholson 
method. The gravitational potential, in turn, is described by Poisson's equa- 
tion, that can be solved using the relaxation method. Our code combines 
these two methods to study the time evolution of a self-gravitating BEC. The 
inefficiency of the relaxation method is balanced by the fact that in subse- 
quent time iterations, previously computed values of the gravitational field 
serve as very good initial estimates. The code is robust (as evidenced by its 
stability on coarse grids) and efficient enough to simulate the evolution of a 
system over the course of 10 9 years using a finer (100 x 100 x 100) spatial grid, 
in less than a day of processor time on a contemporary desktop computer. 

Keywords: 
1. Introduction 

The rotation of spiral galaxies does not follow simple predictions based on 
Kepler's laws. Instead, the rotational velocity curve of most spiral galaxies, 
plotted as a function of radial distance from the galaxy center, remains "flat" 
for a broad range of radii. The standard proposal to resolve this problem is 
to presume the existence of a "dark matter halo", which contains most of 
the mass of a spiral galaxy. To maintain consistency with the predictions of 
the most broadly accepted cosmological models, this halo must necessarily 
consist of "exotic" matter, i.e., matter predominantly composed of something 
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other than baryons. The halo must also be collisionless and not interacting 
with baryonic matter [lj]. 

The existence of such a halo with a suitable geometry can account for 
the observed rotation curves of visible matter. However, a difficult problem 
is to construct a dark matter halo that is gravitationally stable and does not 
predict excessive dark matter densities in the inner parts of the galaxy where 
most visible matter resides. This issue is known as the "cuspy halo problem" 
in the relevant literature 2|. 

A recent proposal 0, |i Ja, @, 0] addresses the cusp problem by a dark mat- 
ter halo that forms a Bose-Einstein condensate (BEC) 0, Ha|. The dynamics 
of a BEC halo may be determined by the balance of the attractive force of 
gravity and a repulsive effective long-range interaction [lo|, 11, 12, [l3[ (see 
also |l4j). As the dark matter halo dominates the gravitational field of a spi- 
ral galaxy, a simulation that is restricted to just the halo should be sufficient 
to determine if a field can be obtained that yields the desired circular orbital 
velocities. 

In the present paper, we discuss a simulation tool that we constructed 
to explore the dynamics of a galactic BEC halo. Our work is based primar- 
ily on our previous simulation of BEC in laboratory conditions, described 
by the non-linear Schrodinger equation, also known in the literature as the 
Gross-Pitaevskii equation. Whereas in the laboratory, a BEC characterized 
by a repulsive interaction is held together by an artificially introduced trap- 
ping potential, in the case of a galaxy floating in empty space, the trapping 
potential must be replaced by self-gravity. A numerical solution must, there- 
fore, simultaneously address the initial value problem of the Gross-Pitaevskii 
equation and the boundary condition problem of Poisson's equation. 

In Sec. [U we introduce the dimensionless form of the Gross-Pitaevskii 
equation used in our computations, and the method used to solve this equa- 
tion efficiently. In Sec. [3] we discuss Poisson's equation for gravity and the 
relaxation method. In Sec. H] we elaborate on the use of physical units that 
are suitable for such a simulation in an astrophysical context. The problem 
of using suitable initial conditions to form a stable halo is briefly discussed in 
Sec. In Sec. El we discuss the implementation of our method in FORTRAN 
and C++, and also comment on the possible use of GPUs for accelerated 
computation. The actual programs are summarized in Sec. [7J Finally, our 
conclusions and outlook are presented in Sec. El 
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2. Solving the Gross-Pitaevskii equation 



A self-interacting, optionally rotating Bose-Einstein condensate is de- 
scribed accurately by a form of the time-dependent nonlinear Schrodinger 



equation known as the Gross-Pitaevskii equation [15|, [16| . For computational 



purposes, it is advantageous to use a dimensionless form of this equation, 



which takes the form 17 



(i-7) 



dip 
~dt 



Hip, 



(1) 



where 7 is a softening parameter (7 = is a valid choice), ip is the wave 
function, t is time, and H is the Hamilton-operator, which in turn is given 
by 

H = -iv 2 + V. (2) 



The potential V is the sum of classical potentials (e.g., gravitational poten- 
tial, trapping potential), the chemical potential, the non-linear term, and a 
rotational term: 

V = <f> + n\ip\ 2 - fi - VL Z , (3) 



where k represents the interaction strength, Q is the angular velocity, and 
L z = i(xd y — yd x ). We assume that the condensate's net rotation is in the 
x — y plane. 



In earlier work 1181, 19, 20, 21 



we solved the Gross-Pitaevskii equation 
can be solved numerically using the Crank-Nicholson method in combination 
with Cayley's formula [22|, in the presence of an isotropic trapping pot ential 
(for a numerical solution in the presence of an anisotropic trap, see 23l. |24| ) . 
In particular, the use of Cayley's formula ensures that the numerical solution 
remains stable and the unitarity of the wavefunction is maintained. 

The value ip n +i of the wavefunction at the (n+l)-th time step is obtained 
from the known values ip n at the n-th time step by solving the following 
equation: 



1 + -iAtH ) ip 



'71+1 



1 



(4) 



After evaluating the right-hand side given ip n , the left-hand side can be solved 
for. If if is a linear operator, this is a linear system of equations for the 
unknown values ip n +i- 
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In the one-dimensional case, the Hamilton operator reads 



H 



1&_ 

2 dx' 



+ V. 



The second derivative can be approximated as a finite difference: 

d 2 tp _ ipk-i ~ 2^fc + tpk+i 
dx 2 (Ax) 2 

Substituting this into Eq. (TjJ, we obtain 



(5) 



(6) 



1 - 



iAt 



V + 



(Ax) 2 
iAt 



iAt 
4(Ax) 



te 1 + tt 1 ) 



(Ax) 



iAt 
4(Ax) 



fc n -i + ^ + i)- (7) 



This is a linear system of equations for the values of ip n+1 on the left-hand 
side, if the values of ip n on the right-hand side are known. Moreover, the 
form of this system of equations is tridiagonal, which can be solved highly 



efficiently using the Thomas algorithm [22 



In the three-dimensional case, one could proceed with solving for ip n+1 
directly, but as the system of equations is no longer tridiagonal, the efficiency 
related to tridiagonal systems is lost. This is why it is preferable to use the 
alternating- direction implicit method, calculating the one-dimensional solu- 
tion in the x, y and z directions, using \At for the time step and \V for the 
potential. This approach is possible because the Hamilton operator can be 
viewed as a sum of three operators, H = H x + H y + H z allowing us to solve 
numerically using fractional time steps as follows: 



1 

~2dx 2 + 3 ' 
1 

2dy 2 + 3 ' 
2 dz 2 3 



(8) 



Substituting these into Eq. ^ (that is, replacing V with V/3 and Ax, 
respectively, with Ax, Ay or Az), using a time step of At/3 we obtain the 
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three fractional iteration steps that correspond to a full iteration with time 
step At. 

As to the nonlinear term, it can be dealt with by a simple iteration 
that converges rapidly. Specifically, we calculate the non-linear term on the 
left-hand side by substituting n in place of <p n+1 and solve the system of 
equations; we then use this solution to recalculate the non-linear term and 
solve again, until convergence is obtained. Because in our case, the non-linear 
term is a cubic term, convergence is very rapid. 



3. Solving Poisson's equation 

The (non-relativistic) gravitational field corresponding to a distribution of 
matter characterized by density p is given by Poisson's equation for gravity: 

V 2 = AirGp, (9) 

where G is the gravitational constant. For a BEC, the mass density is given 
by p = \i^\ 2 m, where m is the mass of the BEC particle. When the BEC 
condensate is described using the dimensionless Gross-Pitaevskii equation, 
m — 1. 

A moderately efficient numerical method for solving Poisson's equation 
is the relaxation method. This method is based on the finite differences 
approximation of the second derivative in Poisson's equation: 



V 2 



<9 2 d 2 <f) d 2 



+ 



dx 2 dy 2 dz 2 

(f>(x — Ax, y, z) — 20(x, y, z) + (j>(x + Ax, y, z) 
Ax~ 2 

(f)(x, y - Ay, z) - 2<f)(x, y, z) + (f)(x, y + Ay, z) 



Ay 



2 



+ (l>(x, y,z - Az) - 2(p(x, y, z) + <p(x, y, z + Az) | Q 2 

Az 2 

where A 2 is the largest of Ax 2 , Ay 2 and Az 2 . If the left-hand side of this 
equation is given (e.g., by Poisson's equation), this equation can be solved 
for <f)(x,y,z). The essence of the relaxation method is the realization that 
these solutions provide successively refined approximations of (f>(x,y,z). In 
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other words, we obtain the iteration formula 



(f) k+1 (x,y,z) = 



Ax 2 Ay 2 Az 2 



+ 



+ 



2(Ax 2 Ay 2 + Ay 2 Az 2 + Az 2 Ax 2 ) 
(f) k (x - Ax, y, z) + (j) k (x + Ax, y, z) 
Ax 2 

<f) k (x, y - Ay, z) + <p k (x, y + Ay, z) 
Ay 2 

4> k (x, y,z- Az) + <f) k (x, y,z + Az) 
Az 2 



4nG\ 



• (11) 



This method is accurate but its convergence is slow. However, after the 
initial configuration of is determined and we iterate the system to the next 
timestep using the Gross-Pitaevskii equation, ip and, consequently, p will 
change very little. Therefore, using the values of n at timestep n as the 
initial estimate for at timestep (n + 1), very rapid convergence is often 
obtained after just a few iterations. 



4. Choice of units 

In order to put the code presented in this paper to use in an astrophysical 
context, it is necessary to restore dimensional units. 

Use of the dimensionless form of the Gross-Pitaevskii equation amounts 
to choosing units such that h = 1 and also m—\, where m is the mass of the 
BEC particle. Given units of length [L], mass [M] and time [T], the choice 
of h — 1 amounts to 

([L] m) 2 -([M]kg) _ m 2 -kg 

(Pis) " 10 s • (12) 

Conversely, choosing a specific numerical value for G amounts to making the 
choice 

G ^fl^l - 2 = 6.67 x lO- 11 -^. (13) 
([M] kg) • ([T] s) 2 kg-s 2 1 ; 

After choosing a specific value for [L] , these two equations allow us to deter- 
mine [M] and [T]: 

[T] = 5.3 x 10 u [Lf 3 ^G s, (14) 
[M] = 5.3 x W^m^VG kg. (15) 
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Specifically, if we choose our unit of length to be 1 kpc ~3x 10 19 m, we get 

[T] = 1.5 x 10 47 ^G s, (16) 

[M] = 1.7 x KT 26 ^ kg ~ 10^G GeV. (17) 

Choosing G = 10~ 100 yields the units 

[T] = 7 x 10 13 s ~ 2.2 x 10 6 year, (18) 

[M] = 8 x 10" 60 kg, (19) 

corresponding to a BEC particle mass of 4.4 x 10 -24 eV. 

Finally, velocity is measured in units of kpc/(2.2 million years) ~ 430 km/s. 



5. Initial and boundary conditions 

Numerically solving a system of coupled differential equations requires a 
set of initial and boundary conditions. 

Specifically, solving the Gross-Pitaevskii equation requires an initial field 
configuration ip(x,y,z) at t = 0. In turn, the numerical solution of Pois- 
son's equation for gravity needs boundary conditions in the form of val- 
ues (f)(x min ,y,z), (f)(x maK ,y,z), <j)(x, y min , z), cj)(x, y max , z), 4>(x,y,z min ) and 
<p(x,y, z max ). 

For the purpose of testing our code, we chose a very simple initial density 
profile, related to solutions of the Lane-Emden equation, centered around 
the origin at x = 0, y = 0, z = 0: 

p = p (R 2 -r 2 ), (20) 

where r 2 = x 2 + y 2 + z 2 and R is a characteristic radius. Given p, we calculate 
ip = y/p, making the initial value of ip purely real everywhere. 

We emphasize that this choice does not necessarily represent a physically 
viable configuration; it was strictly used for code testing and validation. 
Application of the code involves, among other things, choosing initial density 
profiles that reflect valid physical assumptions. For instance, one may opt to 



use a Navarro- Frank- White distribution [25J to model the initial halo density 
at t = 0. 

As to the boundary condition, we simply assume that the gravitational 
field vanishes on the boundaries of the simulation volume. Since in the 
context of Newtonian gravity, the gravitational field is indetermine up to an 
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additive constant, this amounts to the assumption the gravitational field is 
constant on the boundary. While this is somewhat artificial, one may justify 
this choice by noting that real galaxies exist in an external gravitational field 
that in turn is determined by other, more distant galaxies and clusters; this 
external field is imposed upon the dynamics of a real, physical galaxy the 
same way we impose our boundary condition on the model galaxy In any 
case, if the simulation volume is sufficiently large compared to the galaxy 
being simulated, the geometry of the boundary will play no significant role 
in the model galaxy's evolution. 



6. Implementation and testing 

Our code was originally implemented in FORTRAN 90. Later, however, 
we decided to port the code to C++, which resulted in a twofold performance 
improvement. 

We also experimented with a GPGPU0 accelerated version that was de- 
signed to run using the OPENCL library, on an AMD 6900-class graphics 
card. This version yielded considerable improvement in performance, but 
only when single-precision arithmetic was used. The magnitude of the quan- 
tities needed in the astrophysical context (e.g., G = 10 -100 ) precluded the 
use of single-precision floats. The performance improvement of the GPGPU 
version was not sufficient to justify further rewriting the code; therefore, the 
GPGPU implementation was, for the time being, abandoned, although the 
code remains functional. 

The C++ version of the code was heavily tested using various grid sizes. 
When testing the code, we took to heart the important advice offered by the 



authors of Ref. [22J: "You should always first run your programs on very 
small grids, e.g., 8x8, even though the resulting accuracy is [...] poor [...] 
[N]ew instabilities sometimes do show up on larger grids, but old instabilities 
never (in our experience) just go away." We determined that our code runs 
very well on a grid size of 20 x 20 x 20, and the simulation is very rapid; this 
makes it easy to select physically interesting test cases that can be further 
analyzed using a much finer grid. On the other hand, we found that even 
with a grid size of 120 x 120 x 120, a simulation that models the evolution 
over 10 9 years can be completed in the course of a day or so on modern 
desktop computer hardware. 



General Purpose Graphics Processing Unit 
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We specifically tested the numerical robustness of our code by verify- 
ing that the wavefunction if) remains unitary. Even after 50,000 iterations, 
J v \ip\ 2 dV, when evaluated for the entire simulation volume, remained within 
a few percent of its initial value. 

Specifically, we ran a test simulation in a 100 x 100 x 100 kpc 3 volume, 
using an 80 x 80 x 80 cells, using a time step of At = 0.01 units of time 
(corresponding to ~ 22,000 years.) For the entire volume, J \ip\ 2 dV = 
1.7 x 10 102 , corresponding to a total BEC mass of ~ 6.8 x 10 12 M . The 
total simulation time of 50,000 iterations corresponds to ~ 1.1 x 10 9 years; 
the results of this simulation are shown in Fig. [TJ 

7. Program summary 

Our self-gravitating BEC code has several components. A stand-alone 
executable compiled from FORTRAN (bec3p.f90) or C++ (bec3p.c++) 
source performs the simulation. User-definable parameters are in header files 
(parameters3 . f 90 or parameters3 .h). A change in parameters necessitates 
recompiling the executables, but this is quickly accomplished using the sup- 
plied make file (Makefile). The simulation code creates data files at set 
intervals, specifically formatted to facilitate plotting using the gnuplot util- 
ity. A UNIX style shell script (animate . sh) is provided that uses gnuplot 
to create animations in the form of animated GIF files. Finally, yet another 
shell script (video . sh) can be used to convert these GIF files into MPEG-4 
video animations. 

The structure of the main source file (bec3p.f90 or bec3p.c++) is as 
follows: 

• The main program (or _tmain function in the C++ imlementation) 
initializes parameters and executes the main simulation loop. 

• The function init creates the initial distribution. This is the only user- 
modifiable subroutine in the main program file; in the test implementa- 
tion, it creates a naive initial distribution inspired by the Lane-Emden 
equation. 

• The function get_U computes the potential term (without rotation) in 
the Gross-Pitaevskii equation, as given in Eq. 

• The functions calc_rhs_x, calc_rhs_y and calc_rhs_z compute the 
right-hand side of the Gross-Pitaevskii equation (JT]). 
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• The function solve_x, solve_y and solve_z solve for the left-hand side 
of the Gross-Pitaevskii equation according to the scheme represented 
by Eq. (gj. 

• The function thomas solves a tridiagonal system of linear equations 
using the Thomas algorithm. 

• The function getjiormsimp obtains a numerical approximation of J v \ip\ 2 dV 
using Simpson's formula. 

• The functions movieX, movieY and movieZ generate output after a set 
number of iterations (controlled by the variable nstep2) that is set in 
parameters3.h or parameters3 . f 90 in a format that is suitable for 
use with gnuplot. 

• The function get_phi computes the gravitational potential correspond- 
ing to the current values of ip using the relaxation method. 

• The function get_density computes the density \ip\ 2 - 
Output files produced by the program include the following: 

• The files densZnnnnnnn.dat contain values of density (IV'I 2 ), sorted 
and grouped to facilitate plotting with gnuplot in the plane perpendic- 
ular to the Z-axis (i.e., the XY plane). Similarly, the files densXnnnnnnn.dat 
and densY nnnnnnn . dat contain density values sorted and grouped for 
plotting in the plane perpendicular to the X and y-axis. 

• The files gravZnnnnnnn . dat, gravXnnnnnnn . dat and gravYnnnnnnn . dat 
contain values of the gravitational potential, arranged similarly. 

• The files phasZnnnnnnn . dat, phasXnnnnnnn . dat and phasYnnnnnnn . dat 
contain values of the complex phase of ip, again arranged as before. 

The shell script animate, sh provides simple visualization. By default, 
this script create an animated GIF of density values, in the plane perpendic- 
ular to the Z-axis and across the origin, with the color range normalized to 
the 1000th iteration step. Running animate. sh -h prints a brief help message 
that lists all the options of this script; in particular, the option -p that allows 
the user to change the plane for visualization and the option -t that allows 
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the user to change the plot type. The script can also output an animated plot 
of the estimated Keplerian rotational velocity (option -t vrot), for circular 
orbits centered around the origin in the XY plane. 

Finally, the shell script video . sh can be used to convert animated GIF 
output into MPEG-4 video. For this script to work, the utilities f f mpeg and 
mplayer must be installed on the workstation where the script is being used. 

8. Conclusions 

Development of the software code described in this paper is now complete, 
with the code yielding expected results for test cases. The next step is to 
find suitable initial conditions to model a BEC dark matter halo that may 
surround a real galaxy. The question is whether halo configurations can be 
found that are stable over time scales of 10 10 years, and yield circular orbital 
velocities that remain approximately constant at different radii. 

Results of this on-going investigation will be reported when they become 
available. 
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Figure 1: Illustrative example of the evolution of a rotating, self-gravitating BEC, in a 
100 x 100 x 100 kpc 3 volume, simulated using an 80 x 80 x 80 spatial grid. The calculated 
rotational velocity and density are shown after the 1, 250 and 500 iteration time units, the 
latter corresponding to ~ 1 Gyr. The total mass of this condensate is ~ 6.8 x 10 12 Mq. 
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